%% Wage schedule on the grid    
    
    % log wage grid
    nq_cum = cumsum(eqm.nq)-eqm.nq/2; %depends on eqm.nq   
    lnwage_grid = icdf('normal', nq_cum, param.wage_mu(1), param.wage_sigma(1)); %inverse cdf 
    
    % adjust extreme value
    lnwage_grid(eqm.first_nonzero) = lnwage_grid(eqm.first_nonzero+1) + (lnwage_grid(eqm.first_nonzero+1) - lnwage_grid(eqm.first_nonzero+2));   
    lnwage_grid(eqm.last_nonzero) = lnwage_grid(eqm.last_nonzero-1) + (lnwage_grid(eqm.last_nonzero-1) - lnwage_grid(eqm.last_nonzero-2));
    
    % flat at top and bottom
    lnwage_grid(const.Qgrid<eqm.first_nonzero_q)=lnwage_grid(eqm.first_nonzero);
    lnwage_grid(const.Qgrid>eqm.last_nonzero_q)=lnwage_grid(eqm.last_nonzero);
    
    % save
    const.lnwage_grid = lnwage_grid;
    const.wage_grid   = exp(lnwage_grid);  
    
    clearvars nq_cum lnwage_grid

    
%% Wage schedule off the grid

    % Approximate w(q)
    pp_lnwage = interp1(log(const.Qgrid(eqm.nq~=0)), log(const.wage_grid(eqm.nq~=0)), 'spline', 'pp');
    hat_lnwage = @(lnq) ppval(pp_lnwage, lnq);
    hat_wage = @(q) exp(hat_lnwage(log(q)));   
    
    % adjustment: flat at top & bottom
    hat_wage_adj = @(q) (q<eqm.first_nonzero_q).*hat_wage(eqm.first_nonzero_q) ...
                      + (eqm.first_nonzero_q<=q & q<=eqm.last_nonzero_q).*hat_wage(q) ...
                      + (eqm.last_nonzero_q<q).*hat_wage(eqm.last_nonzero_q);
    
    % save
    const.wage_schedule = hat_wage_adj;                              
    const.nq_label_w    = hat_wage_adj([1:8]); %n(q) graph label
    const.nq_label_lnw  = log(hat_wage_adj([1:8])); %n(q) graph label
    
    clearvars pp_lnwage hat_lnwage hat_wage hat_wage_adj

    